function SAT_SST_Step_01_pair_with_HadSST4(P)
    
    % *********************************************************************
    % Load GHCNmV4
    % *********************************************************************
    Data = SAT_SST_func_load_station_data(P);

    % *********************************************************************
    % Pair HadSST4
    % *********************************************************************
    Data.HadSST4_anm  = SAT_SST_func_pair_HadSST4(Data.lon,Data.lat,P,0);

    % *********************************************************************
    % Save data
    % *********************************************************************
    Data.T            = permute(Data.T,          [3 1 2]);
    Data.HadSST4_anm  = permute(Data.HadSST4_anm,[3 1 2]);

    % Subset stations that have HadSST4 paired ----------------------------
    %                         and at leSAT 120 months of observational data
    l_use_1 = ~all(isnan(Data.HadSST4_anm(:,:)),2) & (nansum(~isnan(Data.T(:,:)),2)>=120);
    
    % Subset only high quality stations in further analysis
    l_use_2 = SAT_SST_func_get_high_quality_stations(Data.T,Data.HadSST4_anm,P);
    Data    = struct_subset(Data,l_use_1 & l_use_2,1);
    clear('l_use_1','l_use_2')

    % save data -----------------------------------------------------------
    raw_SAT_full  = Data.T;
    raw_SST_anm   = Data.HadSST4_anm;
    lon           = Data.lon;
    lat           = Data.lat;
    stations      = Data.stations;
    dis           = Data.dis;
    land_fraction = Data.land_fraction;
    file_save = [SAT_SST_func_IO('data',P),'Raw_SAT_SST_',P.date,'.mat'];
    save(file_save,'raw_SAT_full','raw_SST_anm','lon','lat','stations',...
                       'dis','land_fraction','Data','-v7.3');
end
